package br.edu.ufcg.msnlab2.grupo03.romberg;

import br.edu.ufcg.msnlab2.MSNLabException;
import br.edu.ufcg.msnlab2.misc.Function;

/**
 * Implements the Romberg's method of integration
 * 
 * To use this class, just create an object representing an integration problem and
 * call 'run' method.
 * 
 * @author Thiago Santos
 * @author Joao Paulo
 */
public class RombergMethod {

	Function f;
	double lowLimit;
	double upLimit;
	double[] previous;
	double[] current;

	/**
	 * Creates a new RombergMethod object.
	 * 
	 * All parameters are required and lowLimit cannot be greater than upLimit
	 * 
	 * @param f the function to be integrated
	 * @param lowLimit the lower limit of the integration
	 * @param upLimit the upper limit of the integration
	 * @throws MSNLabException 
	 */
	public RombergMethod(Function f, double lowLimit, double upLimit) throws MSNLabException {
		if (f == null) {
			throw new MSNLabException("A function must be specified");
		}
		if (lowLimit > upLimit) {
			throw new MSNLabException("The inferior limit cannot be greater than the superior limit");
		}
		this.f = f;
		this.lowLimit = lowLimit;
		this.upLimit = upLimit;
	}
	
	public RombergMethod(){
		
	}

	private void resetData() {
		previous = new double[1];
		current = null;
	}

	/**
	 * Applies the trapezium rule for a function, from lowLimit to upLimit.
	 * 
	 * @param f function to apply the rule to
	 * @param lowLimit inferior limit
	 * @param upLimit superior limit
	 * @return
	 */
	public static double trapeziumRule (Function f, double lowLimit, double upLimit) {
		return ((upLimit - lowLimit) * (f.evaluate(lowLimit) + f.evaluate(upLimit)))/2.0;
	}

	/**
	 * Static helper function that instantiates a RombergMethod object and runs Romberg's method.
	 * @throws MSNLabException 
	 * @see constructor
	 * @see run
	 */
	public static double execute (Function f, double lowLimit, double upLimit, int maxIterations, double tolerance) throws MSNLabException {
		return new RombergMethod(f, lowLimit, upLimit).run(maxIterations, tolerance);
	}

	/**
	 * Executes the Romberg's Integration method on this object with 
	 * a limited number of iterations or until the tolerance is 
	 * achieved.
	 * 
	 * The tolerance is defined by a threshold between two subsequent approximations,
	 * if the difference between them is greater then the tolerance, the next approximation
	 * is calculated. It is stopped otherwise.
	 * 
	 * @param maxIterations the maximum number of iterations to perform
	 * @param tolerance the minimum difference between two subsequent approximations
	 * @return the integration result
	 * @throws MSNLabException 
	 */
	public double run (int maxIterations, double tolerance) throws MSNLabException {
		if (tolerance < 0)
			throw new MSNLabException("Tolerance must be non-negative");
		resetData();
		if (lowLimit == upLimit) return 0;
		if (maxIterations == 0) return 0;
		previous[0] = trapeziumRule (f, lowLimit, upLimit);
		int iteration = 0;
		while (true) {
			iteration++;
			current = new double[iteration+1];
			if (iteration > maxIterations) {
				break;
			}

			/* define the first value of this line based on the first value
			 * of the last one */
			current[0] = romberg1(iteration);

			for (int m = 1; m < current.length; m++) {
				current[m] = romberg2(m);
			}

			previous = current;
			current = null;
			
			/* check tolerance */
			if (Math.abs(previous[previous.length-1] - previous[previous.length-2]) <= tolerance) {
				break;
			}
		}
		return previous[previous.length - 1];
	}
	
	static double heightN(double lowLimit, double upLimit, int n){
		return (upLimit - lowLimit)/Math.pow(2, n);
	}

	double romberg1(int iteration){
		double sum = 0;
		for (int i = 1; i <= Math.pow(2, iteration-1); i++){
			double x = lowLimit + ((2*i) - 1)*heightN(lowLimit, upLimit, iteration);
			sum += f.evaluate(x);
		}
		return (previous[0])/2 + heightN(lowLimit, upLimit, iteration)*sum;
	}
	
	double romberg2(int m) {
		double r1 = current[m-1];
		double r2 = previous[m-1];
		double pow = Math.pow(4, m);
		return r1 + (r1 - r2)/(pow- 1);
	}
	
}
